﻿#define _CRT_SECURE_NO_WARNINGS
#pragma execution_character_set("utf-8")

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <string.h>
#include <windows.h>  // 新增头文件


#define MAX_DATA_SIZE 200// 最大数据行数
#define MAX_ITER 1000  // Levenberg-Marquardt算法的最大迭代次数
#define EPSILON 1e-9   // 收敛判定阈值（误差变化小于此值时停止迭代）

// 数据结构体：存储单个受试者的试验数据
typedef struct {
    double time;// 给药时间（小时，h）
    double conc_plasma;// 血药浓度观测值（ng/mL）
    double skin_thickness;// 皮肤厚度（cm，从原始数据mm转换而来）
    double BMI;// 身体质量指数（kg/m²）
} DataPoint;

/**
 * 解析带有“±”符号的数值（如"0.02��0.0"），提取均值
 * @param str 输入字符串
 * @return 解析后的浮点数值（如0.02）
 */
double parse_value(const char* str) {
    char buffer[32];
    strncpy(buffer, str, sizeof(buffer) - 1);
    buffer[sizeof(buffer) - 1] = '\0';
    char* sep = strstr(buffer, "��");
    if (sep) *sep = '\0';
    return atof(buffer);
}

/**
 * 安全读取CSV文件，解析字段并填充到DataPoint数组
 * @param filename CSV文件名
 * @param data 存储数据的数组
 * @return 成功读取的数据行数，失败返回-1
 */
int read_csv(const char* filename, DataPoint* data) {
    FILE* file;
    errno_t err = fopen_s(&file, filename, "r");
    if (err != 0 || file == NULL) {
        perror("无法打开文件");
        return -1;
    }

    char line[512];
    fgets(line, sizeof(line), file);// 跳过标题行（首行）

    int count = 0;
    while (fgets(line, sizeof(line), file) && count < MAX_DATA_SIZE) {
        // 定义临时变量存储各字段
        char id[10], age[10], bmi[10], skin_thick[10], time[10], plasma[20], surf[20], deep[20], effect[10];
        // 使用sscanf_s安全解析CSV行，指定每个字段的缓冲区大小
        if (sscanf_s(line, "%[^,],%[^,],%[^,],%[^,],%[^,],%[^,],%[^,],%[^,],%s",
            id, sizeof(id), age, sizeof(age), bmi, sizeof(bmi), skin_thick, sizeof(skin_thick),
            time, sizeof(time), plasma, sizeof(plasma), surf, sizeof(surf), deep, sizeof(deep),
            effect, sizeof(effect)) != 9) {
            continue; // 解析失败时跳过该行
        }

        data[count].time = atof(time); // 给药时间
        data[count].conc_plasma = parse_value(plasma);// 血药浓度（处理“��”）
        data[count].skin_thickness = atof(skin_thick) / 10.0;// 皮肤厚度（mm→cm）
        data[count].BMI = atof(bmi);// BMI
        count++;
    }

    fclose(file);
    return count;// 返回有效数据行数
}

/**
 * 药物动力学模型：描述药物在皮肤层扩散和血药浓度代谢
 * @param D 扩散系数（cm²/h）
 * @param k_abs 吸收速率常数（1/h）
 * @param k_elim 消除速率常数（1/h）
 * @param Vd 分布容积（mL）
 * @param t 给药时间（h）
 * @param L 皮肤厚度（cm）
 * @return 预测的血药浓度（ng/mL）
 */
double model(double D, double k_abs, double k_elim, double Vd, double t, double L) {
    double c_skin = D * L * (1 - exp(-k_abs * t));// 皮肤层药物浓度（μg/cm²）
    return (c_skin * k_abs / Vd) * exp(-k_elim * t);// 血药浓度（ng/mL）
}

/**
 * Levenberg-Marquardt算法优化模型参数
 * @param params 待优化参数数组：[D, k_abs, k_elim, Vd]
 * @param data 试验数据数组
 * @param n 数据量
 */
void lm_optimize(double* params, DataPoint* data, int n) {
    double lambda = 0.01;// 阻尼因子初始值
    double J[MAX_DATA_SIZE][4], r[MAX_DATA_SIZE]; // 雅可比矩阵（n行×4列，// 残差数组
    double prev_error = INFINITY;// 前一次迭代的误差

    for (int iter = 0; iter < MAX_ITER; iter++) {
        double error = 0.0;
        // 计算残差和雅可比矩阵
        for (int i = 0; i < n; i++) {
            double pred = model(params[0], params[1], params[2], params[3],
                data[i].time, data[i].skin_thickness);// 模型预测值
            r[i] = pred - data[i].conc_plasma;// 残差（预测值-观测值）
            error += r[i] * r[i]; // 总误差（平方和）

            // 数值法计算雅可比矩阵（偏导数）
            double h = 1e-6;
            for (int j = 0; j < 4; j++) {
                double tmp = params[j];// 保存原参数值
                params[j] += h;// 微小扰动
                J[i][j] = (model(params[0], params[1], params[2], params[3],
                    data[i].time, data[i].skin_thickness) - pred) / h;// 计算偏导
                params[j] = tmp; // 恢复参数
            }
        }

        // 完整矩阵运算， 构建正规方程 JTJ 和 JTr
        double JTJ[4][4] = { 0 }, JTr[4] = { 0 };
        for (int j = 0; j < 4; j++) {
            for (int k = 0; k < 4; k++) {
                for (int i = 0; i < n; i++) {
                    JTJ[j][k] += J[i][j] * J[i][k];
                }
            }
            JTJ[j][j] += lambda;
            for (int i = 0; i < n; i++) {
                JTr[j] += J[i][j] * r[i];
            }
        }

        // 解方程（简化版）
        double delta[4] = { 0 };
        for (int j = 0; j < 4; j++) {
            delta[j] = -JTr[j] / (JTJ[j][j] + 1e-6);// 参数增量
            params[j] += delta[j];// 更新参数
            if (params[j] < 0) params[j] = 0;// 强制非负约束
        }

        // 动态调整 lambda
        if (error < prev_error) {   // 误差减小，减小lambda（更接近高斯-牛顿法）
            lambda *= 0.1;
            prev_error = error;
        }
        else {
            lambda *= 10.0;  // 误差增大，增大lambda（更接近梯度下降法）
        }

        if (fabs(prev_error - error) < EPSILON) break;
    }
}

int main() {
    DataPoint data[MAX_DATA_SIZE];// 存储所有试验数据
    int n = read_csv("clinical_trial_data.csv", data);// 读取CSV文件
    if (n <= 0) return 1;

    // 初始参数猜测：扩散系数、吸收速率、消除速率、分布容积
    double params[4] = { 1e-6, 0.01, 0.01, 1000.0 };  // 更合理的初始猜测

    // 执行Levenberg-Marquardt优化
    lm_optimize(params, data, n);

    SetConsoleOutputCP(65001);  // 设置控制台输出编码为UTF-8

    printf("优化后的参数:\n");
    printf("扩散系数D = %.4e cm²/h\n", params[0]);
    printf("吸收速率k_abs = %.4f 1/h\n", params[1]);
    printf("消除速率k_elim = %.4f 1/h\n", params[2]);
    printf("分布容积Vd = %.2f mL\n", params[3]);

    return 0;
}